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Abstract 



We study a one-dimensional fixed-energy version (that is, with no input 
or loss of particles), of Manna's stochastic sandpile model. The system has 
a continuous transition to an absorbing state at a critical value Q c of the 
particle density. Critical exponents are obtained from extensive simulations, 
which treat both stationary and transient properties. In contrast with other 
one-dimensional sandpiles, the model appears to exhibit finite-size scaling, 
though anomalies exist in the scaling of relaxation times and in the approach 
to the stationary state. The latter appear to depend strongly on the nature 
of the initial configuration. The critical exponents differ from those expected 
at a linear interface depinning transition in a medium with point disorder, 
and from those of directed percolation. 
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I. INTRODUCTION 



Sandpile models are the prime example of self-organized criticality (SOC), or scale- 
invariance in the apparent absence of tuning parameters SOC in a slowly driven 

sandpile has an analog in the absorbing-state phase transition of the corresponding nondriven 
or fixed-energy sandpile (FES) While most studies of sandpiles have probed the slow- 

driving limit (addition and loss of sand grains at an infinitesimal rate), there is great interest 
in understanding the scaling properties of FES models as well [J^-jlT]] . One may hope that 
these systems, which are isolated and translation-invariant, will prove easier to characterize 
than driven sandpiles, which have eluded a rather massive theoretical and computational 
effort to obtain precise critical exponents or assign universality classes. In this paper we 
present extensive numerical results for a particularly simple one-dimensional FES. Analogous 
studies of higher dimensional systems can be found in |7||3|,|i"2"| (A review of these results and 
related models can be found in Refs. 0,0). 

One of the central features of sandpile models is the presence of a conserved field, which 
may be thought of as the density of sand grains, or, in our nomenclature, the energy density. 
This field is conserved locally as well as globally, and couples to the activity density, which 
is the order parameter. The configuration evolves through a series of "toppling" events, i.e., 
redistributions of energy from an active site to its neighbors, which may be either determinis- 
tic or stochastic. The well known Bak-Tang-Wiesenfeld (BTW) sandpile has a deterministic 
toppling rule, and is also Abelian (the final absorbing configuration is independent of the 
order of toppling), which allows many stationary-state properties of the driven sandpile to 
be found exactly |13| , [14| . A less desirable aspect of the deterministic dynamics, on the other 
hand, is that in the steady-state only a small subset of the possible configurations (deter- 
mined by the initial state) are visited fl3f . This leads to strong nonergodic effects in the 



FES version of the BTW automaton |T2[ . Here we study a one-dimensional stochastic FES 
sandpile to elucidate the consequences of stochastic rules. 

In our model, a variant of the Manna sandpile |i5| , |I6|| , the pair of particles liberated when 
a site topples move independently, and with equal probability, to the nearest neighbors on 
the left or right. We may therefore think of the particles as independent random walkers. 
(There is no restriction on the number of walkers at a given site.) While the walkers are 
independent as far as their hopping direction is concerned, their mobility does involve an 
interaction: isolated particles cannot move, but if two or more particles occupy the same 
site, the site may topple. Thus for densities ( = N/L < 1 (N is the number of walkers on 
a ring of L sites), absorbing configurations exist, in which all walkers are paralyzed. But 
since the fraction of absorbing configurations vanishes as £ — > 1, it is reasonable to expect 
a phase transition from an absorbing to an active phase at some ( c < 1. Simulations bear 
this out, and show that there is a continuous transition at ( c ~ 0.9488. 

Our goal in the present paper is to determine the critical behavior of the model. To this 
end, we analyze the static behavior and the approach to the steady state in the language 
of absorbing-state phase transitions, studying the order parameter, its fluctuations, and 
temporal correlations. We also investigate the cumulative local activity, whose dynamics 
obeys a linear diffusion equation in the presence of an effective quenched noise term that 
reflects the stochastic redistribution rules, the Linear Interface Model (LIM) p]-|2"0"||. Our 
main finding is that the one- dimensional Manna model defines a universality class different 
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both from that of directed percolation, and of the linear interface model. The existence of 
long-range correlations that arise from a conservation law is believed to be at the basis of 
this new universality class |f7||3|, |i"2"| , |2~T . 



The paper is structured as follows. In Sec. II we define the model and our simulation al- 
gorithm. Numerical results are analyzed, in the contexts of absorbing-state phase transitions 
and of driven interfaces, in Sec. III. In Sec. IV we summarize and discuss our findings. 



II. MODEL 

The model is defined on a one-dimensional lattice with periodic boundaries. The con- 
figuration is specified by the energy (or number of walkers) Z{ = 0, 1, 2, ... at each site; sites 
with Zi > 2 are said to be active. A Markovian dynamics is defined by the toppling rate, 
which is unity for all active sites, and zero for sites with < 2. When a site i topples, 
Zi — > Zi — 2 and the two particles liberated move independently to randomly chosen nearest 
neighbors j and f (j, f G {i + l,i — 1}). (Thus j = f with probability 1/2.) In most cases 
we use random sequential dynamics: the next site to topple is chosen at random from a list 
of active sites, which must be updated following each toppling event. The time increment 
associated with each toppling is At = 1/Na, where is the number of active sites just 
prior to the toppling. In this way, (N^) sites on average topple per unit time, just as they 
would in a simultaneously updated version of the model, in which all active sites at time 
t n are toppled simultaneously and At = 1. We expect the two dynamics to be equivalent 
insofar as scaling properties are concerned, and note that the latter choice was used in some 
of the interface representation studies discussed below. 

In most of our simulations, the initial condition is generated by distributing (L particles 
randomly among the L sites, yielding an initial (product) distribution that is spatially 
homogeneous and uncorrelated. Once the particles have been placed, the dynamics begins 
(we verified that allowing some toppling events during the insertion phase has no effect on 
the stationary properties). The condition of having at least one active site in the initial 
configuration is trivially satisfied on large lattices, for the ( values of interest, i.e., close to 
the critical value. In fact, for large L, the initial height at a given site is essentially a Poisson 
random variable, and the probability of having no active sites decreases exponentially with 
the lattice size. It is worth remarking that while the initial conditions are statistically 
homogeneous, the energy density is not perfectly smooth. For 1 C I < L, the energy 
density on a set of I sites is essentially a Gaussian random variable with mean ( and 
variance ~ l~ d . Such initial fluctuations relax on a time scale ~ l z ° with zd > 2. That 
is, the exponent governing diffusive relaxation must be at least as large as for independent 
random walkers. (Recall that in our model regions of lower density will exhibit less activity, 
hence a slower redistribution of particles.) The initial value of the critical-site density, 
p c (the density of sites with = 1), moreover, is generally far from its stationary value, 
complicating relaxation to the steady state. These observations are pertinent to the studies 
of modified initial conditions reported in Sec. Hid. 

If after some time the system falls into a configuration with no active sites, the dynamics 
is permanently frozen, i.e., the system has reached an absorbing configuration. We shall 
see that as we vary (, the model exhibits a phase transition separating an absorbing phase 
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(in which the system always falls into an absorbing configuration), from an active phase 
possessing sustained activity. 



III. SIMULATION RESULTS 
A. Absorbing-state phase transition 

The presence of an absorbing-state phase transition is a general feature of FES models 



I"2]j . In analogy with standard analyses of such transitions, it is possible to derive a set 
of Langevin equations describing their critical behavior in high dimensions. Following Ref. 
|12|1 , we can write the dynamical equations for the density of active sites p a and the energy 



density ( as 

d t p a = DV 2 p a + Ppa ~ bp\ + w(p a + p^T] (1) 



dt( = XV 2 p a , (2) 



where /i, b, w and A are coupling constants and r\ is a Gaussian noise. Eqs. (0-@) closely 
resemble the Langevin field theory of directed percolation (DP), apart from an additional 
coupling between the activity and the energy field, which is conserved by the dynamics. Due 
to this coupling, the FES critical behavior is not in the DP universality class, but belongs to 
a new universality class embracing all systems with multiple absorbing states in which the 
order parameter is coupled to a static conserved field [pT| , p2| . (By 'static' we mean that this 
field is frozen in regions where the order parameter is null.) The field theory for a reaction- 
diffusion model in this universality class has been derived exactly; it reproduces in some 
limits Eqs. (|l]-0) ||22|| . It is thus interesting to study the critical behavior in d = 1, to see if it 
is possible to generalize the known results for fixed-energy sandpiles to one-dimensional 
systems. 

We simulated the model using system sizes ranging from L = 100 to about 10 4 sites. 
(Since ( = N/L with N an integer, we are obliged to use different sets of L values to 
study different energy densities (.) In stationary-state simulations, we collect data over an 
interval of t m time units, following a relaxation period of t r . For small system sizes, t m and 
t r are of the order of 10 3 , but for our largest systems (L ~ 10 4 ) we used t r > 5 x 10 6 and 
tm = 2.5 x 10 6 . We verified that our results show no systematic variation with time for 
t > t r . A run consists of N s independent trials, each with a different initial configuration, 
with N s ranging from 2 x 10 5 for L = 100, to 500 or 1000 for L ~ 10 4 . In practice t m is 
limited because for ( ~ ( c , the survival probability decays sensibly over this time scale; in 
some cases only about 25% of the trials survive to time t r + t m . 

Fig. 1 shows the overall dependence of the stationary active-site density as a function 
of C; the points represent extrapolations of results for L = 100 - 5000 to the L — > oo limit. 
The data indicate a continuous transition from an absorbing state (p a = 0) to an active one 
at ( c in the vicinity of 0.95. 

Our first task is to locate the critical value ( c . We studied the stationary active-site 
density, p a , and its second moment p 2 , anticipating that as in other absorbing-state phase 
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transitions, the active-site density (i.e., the order parameter) will obey finite-size scaling 
2|: 

pZ(A,L)=L-^n(L 1 ^A), (3) 

where A = ( — ( c , and 1Z is a scaling function with lZ(x) ~ x 13 for large x, since for 
L ^> £ ~ A~ u± we expect ~ A^ (here £ is the correlation length). When A = we 
have that p~a(0, L) ~ L~ l3 ^ u± . For A > 0, by contrast, approaches a stationary value, 
while for A < it falls off as L- d . Thus in a double-logarithmic plot of p a versus L (see 
Fig. 2), supercritical values (A > 0) are characterized by an upward curvature, while for 
A < the graph curves downward. Using this criterion (specifically, zero curvature in the 
data for L > 1000), we find ( c = 0.94887(7), with the uncertainty reflecting the scatter in 
our numerical results for the curvature (see Fig. 2, inset). The associated exponent ratio 
is P/u± = 0.235(11). A similar analysis of the data for p\ yields ( c = 0.94883(5) with an 
exponent of 2(3 /vx = 0.483(18). We therefore adopt the estimates ( c = 0.94885(7) and 
P/v x = 0.239(11). 

We measured the autocorrelation function for the number of active sites N a , 

c{t) ~ W) - W 2 (4) 

in the stationary state. To obtain clean results for C(t) we study surviving trials in relatively 
long runs (from t m = 2 x 10 5 for L = 625, to 5 x 10 6 for L = 10 4 ; this obliges us to reduce 
our sample to 200 surviving trials for L < 2500 and 100 surviving trials for L > 5000). 
Results for ( = 0.9488 are shown in Fig. 3: C(t) decreases monotonically, but does not 
follow a simple exponential decay. To study the dependence of the relaxation time on 
system size, we determine the temporal rescaling factor r required to obtain a data collapse 
between C(t; L/2) and C(t/r; L). A good collapse is possible (see Fig. 3), but the rescaling 
factor depends on L; for L = 2 n ■ 625, we use t* = t/r n with r = 2.93, 2.91, and 2.89 
for n = 1, 2, and 3, respectively. The rescaling factor r appears to approach a limiting 
(L — > oo) value of 2.80(5), corresponding to a relaxation time that scales as r ~ L v w' Vl - with 
v\\jvx = z = In rj In 2 = 1.49(3). We also studied the half-life Th (the time for the survival 
probability to decrease by a factor of 1/2 in the stationary state) at ( c . This relaxation time 
exhibits a clean power-law dependence on system size, r# ~ L z , with z = 1.42(1). 

Next we examine the stationary scaling of the order parameter away from the critical 
point. We determined the stationary active-site density ~p^((,L), for ( in the vicinity of 
( c = 0.94885, and for system sizes L = 100 - 5000. We analyze these data using the 
finite-size scaling form of Eq. (^), which implies that a plot of L^' Uj -p~^(A, L) versus L 1 ^ U± A 
should exhibit a data collapse. We shift each data set (in a log-log plot of L l3 ^ u± p a versus 
L 1 / U± A), vertically by (fl/vx) InL, using f3/vx = 0.239 as found above, and determine the 
horizontal shift S(L) required for data collapse. These are found to follow the relation 
S(L) = z/J 1 InL, with z/J 1 = 0.553(3). That these values yield an excellent data collapse is 
evident from Fig. 4. The slope of the scaling plot (linear regression using the 25 points with 
\n(L 1 ^ u± A) > —0.5), yields (3 = 0.410(4). This is somewhat smaller than, but consistent 
with, the estimate (3 = 0.43(2) obtained by combining z/J 1 = 0.553(3) and (3/vx = 0.239(11). 
We adopt (3 = 0.42(2) as our final estimate. 
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B. Interface Representation 



Interface representations are useful for studying both sandpiles [|17|-|2"0|| and particle sys- 
tems such as the contact process f30fl . Using such a representation, the activity history of 
the sandpile can be described in terms of an equation of motion that is an exact analog 
of a driven interface advancing in a quenched disordered medium. Our one-dimensional 
case is an example of a larger class of models, for which a mapping can be provided from 
the sandpile to an interface field [|1TH 20[| . Such a description can be useful, for example, in 
understanding the role of microscopic rules and in defining novel quantities for study. 

The mapping is constructed as follows. Let the interface height variable Hi(t) count 
the number of topplings at site i up to time t. The dynamics of Hi(t) follows from the 
observation that the toppling of a site depends (through its height z^) on the number of 
particles it has received, and on the number of times it has toppled. Each time a neighbor 
of site i topples, site i gains one, two, or no particles, with probabilities of 1/2, 1/4 and 1/4, 
respectively. This allows a 'projection trick', in which the dynamics of Hi(t) is constructed 
from particle counting |]20| , p~2|] . Topplings at the sites neighboring i give rise to an average 
incoming flux n™ at site i, plus a fluctuating part dnf 1 , which records, in an exact manner, 
the random choices made for each particle. The average flux minus the number of topplings 
at i gives rise to a discrete Laplacian, while the fluctuating part can be mapped to a noise 
term r(i,H), which can be evaluated for each toppling. Following these prescriptions, we 
find a discrete interface equation with quenched noise, which reproduces the sandpile activity 
2|: 

= 6[V 2 ^ + - 1 + r(i, H)] = e(/ 4 ) (5) 

where G is the step function, so that Hi is nondecreasing. Since Hi(t) is restricted to integer 
values the time derivative should be interpreted as the rate of transitions from H to H + 1. 
The 'force' fi at site i consists of the following contributions. The (discrete) Laplacian term 
arises, as noted above, from the balance between the average number of particles site i gains 
from, and loses to, its neighbors in toppling events. The 'columnar' noise term Fi represents 
the number of particles initiallly at site i, while r(i, H) is the fluctuation part of the particle 
flux associated with the H-th toppling of site i. The step-function character of the interface 
equation is discussed in greater detail in Ref. [^U] . 

The noise-term r(i, H) is a quenched variable that arises from a mapping of the 'annealed' 
disorder of the original sandpile. Note that one could also consider the random choices of 
the particles' destinations at each toppling as a quenched, random landscape that is chosen 
in advance for each sandpile configuration. It is then easier to see that this maps into a 
quenched disorder term for the associated interface equation. This noise term is particle- 
conserving, and it has at each i a random- walk- like correlator ([t(H + AH, i) — t(H, i)] 2 ) ~ 
AH. If these correlations could be neglected, the scaling of the interface would follow that 
of the linear interface model in a quenched random field, as appears to be the case in two- 
and three-dimensional systems |12| . One should note, however, that long range-correlations 
in the noise can change the universality class, as has been shown numerically for another 



interface model pj|]. The possible LIM universality classes are still under debate 132" 



In the interface description the model undergoes a depinning transition at the critical 
force F c = ( c |2|j2|]. At this point power-law correlations develop in the history of top- 
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plings. The interface behavior, assuming simple scaling, is described by two exponents: the 
roughness exponent a, and the early-time exponent fiw Introducing the width W as usual, 

W*(t,L) = ([H l (t)-H(t)} 2 ), (6) 
(here H(t) is the mean height), these exponents are defined via 

w\t,D~[^: ;<;-, ( 7) 

where the crossover time t x ~ L z . Assuming that simple scaling holds and that the cor- 
relations in the interface can be described by a single length scale, we have the exponent 
relation j3w z = ol. 

This scaling picture, familiar from the study of surface growth, was recently shown to 



apply in the case of a simple absorning-state phase transition ]30| . For the one-dimensional 
Manna model the situation is complicated by several factors. First, the noise appearing in 
the interface description has two components, r(i, H) and F^. The interface behavior will 
therefore exhibit a crossover from a regime dominated by the initial configuration (reflected 
in Fj) to a randomness-dominated regime. This effect also appears in higher dimensions, 
but in d — 1, due to the meager phase space, relaxation is much slower and transient effects 
may be much more severe. Note that the if -independent terms on the r.h.s. of the interface 
equation can be understood as an initial height profile, which then relaxes to the asymptotic 
rough state ||20|| . It is easy to see that by varying such an initial condition different transients 
can be generated. 

A problem with one-dimensional models is that the two-point correlation function of the 
surface roughness scales with a different exponent, a; oc , which for fundamental reasons is 
limited to a\ oc < 1. The exponent a can attain a larger value, for example a = 1.25 for the 
one-dimensional linear interface model |26||. These exponents are related via a = a\ oc + k, 



where k measures the divergence of the height difference between neighboring sites with L. 



In the corresponding 1-d linear interface model, this is termed 'anomalous scaling' ||26| . |27 
This means that as t — > oo, the typical height difference between neighboring sites increases 
without limit. Since larger systems have a greater lifetime, this has implications for the 
roughness as measured by W^ at = limt_ +tc W 2 (t, L), with t c the time at which the absorbing 
state with no activity is reached, i.e., the lifetime. The saturation width W sa t scales as L a , 
with a related to j3w and z as above. In models exhibiting an absorbing state, such as the 
contact process or a FES, the width saturates only because all activity eventually ceases; 
the width in surviving trials does not saturate. This is in marked contrast to the behavior 
of interface models, in which the width saturates due to the Laplacian term representing 
surface tension. (Continuum descriptions of absorbing-state phase transitions, and their 
associated interface representations, likewisw include a Laplacian term, but the noise driving 
the interface grows without limit, so saturation is not required while there is activity pCfl .) 

Finally, in absorbing-state models, interface scaling appears to be strongly linked to the 
approach to the stationary state. In a model with simple scaling (i.e., unique diverging 
length and time scales, and no conserved quantities), the growth exponent j3w is related to 
the critical exponent governing the initial decay of activity via f3w + = 1 ||30|| . In the 
present case relaxation is complicated by effects that may mimic (for a certain time) colum- 
nar disorder, and slowly relaxing perturbations (the long- wavelength density fluctuations 
mentioned in Sec. II). 
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In simulations we first studied the time-dependent mean interface width W 2 (t,L) in 
systems of size L = 1253, 2506, 5012, 10024, and 20048, at the critical point ( = 0.94892. 
The dependence of the saturation width on system system size yields a = 1.42(1). We may 
then attempt to collapse the data for W 2 (t, L) using this exponent, and varying z to obtain 
the best collapse; in this way we find z = 1.65(2). The resulting scaling plot (Fig. 5), of 
W 2 = W 2 /L 2a versus t = t/L z , shows a good collapse, and an apparent power-law growth in 
the roughness, following an initial transient. From the scaling relation (3w = a/z we obtain 
Pw — 0.863(13), while a direct fit to the time-dependent width data yields j3w — 0.87(2). 

For comparison, an independent series of studies at £ = 0.9490 were performed to de- 
termine the lifetime t c and saturation width W c for L = 400 to L — 6400. Power-law 
fits to these data yield essentially consistent results, i.e., a = 1.48(2), j3w = 0.86(2) and 
z = 1.70(3). Fig. 6 shows a clear power-law dependence of the saturation width on the 
lifetime in individual runs. This indicates that a normal Family- Vicsek-style scaling plot of 
W(t, L) as in Fig. 5 is problematic, since as noted above the width in any surviving run 
does not saturate as would be the case in normal interface models, and the saturation width 
W ~ w c will therefore depend on how long the particular configuration stays alive. 

Similarly to the case of the linear interface model, we find that there is an independent 
local roughness exponent ai oc that describes the two-point k-th order height-height correla- 
tion function G k {r) = (\H i+r - Hi\ k ) ~ r fca <- for r < f (t) ~ t 1 / 2 . We find a loc = 0.59(3) 
for k = 1 ... 8. This shows that the interface is simply self-affine, not multifractal. We 
note, however, that the interface exponents of the Manna sandpile are not those of the 
one-dimensional LIM. This is most likely due to the fact that, perhaps differently from 
the two- and higher-dimensional cases, here it is important that the noise term increases 
in strength with the propagation of the interface, or with sustained activity. This arises 
from the fact that the fluctuations in the number of grains received at each site, reflected 
in r(H,i), follow random- walk characteristics. Thus the anomaly exponent k indicates an 
even stronger dependence of the step height on L at saturation than in the LIM. The same 
is also true if the step height is considered as a function of time for t < t c . We find that 
^Manna ~ 0.82 ~ Klim + 0.5. Due to this time- or height-dependence of r there is a direct 
coupling of t c and w c with the final average height, H c (see Fig. 7). As discussed above, this 
result follows, in the interface picture, from the fact that the noise strength increases with 
time. In the LIM, t c and w c are not correlated on a sample-to-sample basis. 



C. Initial Relaxation 



A comparison of our results for stationary correlations (z ~ 1.45) with those for the 
interface (z ~ 1.7) suggests that the relaxation to the steady state and relaxation of fluctua- 
tions within the stationary regime show distinct scaling properties. This is supported by our 
results for the relaxation of the active-site density. At the critical point of a simple absorbing- 
state model such as the contact process, the activity density p a exhibits an initial power-law 
decay, p a ~ t~ 9 , followed by a crossover to the quasi-stationary value p a ~ L~PI Ul - P9]| . As 



noted above the growth exponent is related to the activity- decay exponent via 9 + j3w = 1 
if only one timescale is present [ftU|] . A plot of L f3 ^ u± p a (t) versus t/L z yields a data collapse 
to a scaling function that is independent of L. In the present case (Fig. 8), we see that 
the collapse is imperfect, and that the form of p a {t) changes with L. Here z was chosen 
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so as to optimize the collapse at long times, yielding z = 1.75(3). For large systems, the 
active-site density exhibits three distinct regimes before reaching the quasi-stationary state: 
an initial power-law decay (I), followed by a crossover to a slower power-law regime (II), and 
finally a rapid decay (III) to the stationary state. For L = 20048, the exponents associated 
with regimes I and II are 0.163 and 0.144, respectively. While the latter exponent is in 
reasonable agreement with the scaling relation j3w + = 1, it is clear that relaxation to 
the stationary state is more complicated for the sandpile than for, say, the contact process, 
which presents a unique power-law regime. A qualitative explanation may be found in the 
interface representation: the initial dynamics is dominated by relaxation of the initial grain 
profile Zi, which in the interface language means that at short times the columnar noise, F iy 
dominates. 

Thus another facet of the relaxation process is the approach of the mean height to 
its global value, (, at a site with initial height z(0). In Fig. 9 we plot the mean height 
{z(t)\z(0)) in at system of 1400 sites at ( c , averaged over 2000 trials. The inset shows that 
the asymptotic approach to ( is approximately power-law, \(z(t)\z(0)) — CI ~ ~t~^-> with 
=~ 0.46, 0.45, 0.47, 0.50, and 0.53 for z(0) = 0, 1, 2, 3 and 4, respectively. All of these 
exponents are close to = 1/2, the value expected for uncorrelated diffusion. 

D. Effects of initial preparation 

In light of the complicated pattern of relaxation noted in the preceding subsection, it is 
well to recall that initial-decay studies of simple models such as the contact process employ 
a uniform, featureless initial configuration, far from the absorbing state. This is usually 
arranged by setting the the activity density to unity, an option that is not available in the 
present case, if we want to fix the particle density to its critical value ( c . Since the height is a 
discrete variable, some local variations in density and activity are inevitable. All the results 
described up to now were obtained using the random deposition (RD) initial preparation 
described in Sec. II. In this subsection we report on studies using initial configurations in 
which the N = (L particles are distributed so as to reduce density fluctuations. We again 
employ the sequence L = 2 n ■ 1253, N = 2 n - 1189 (n = 0, 4) corresponding to ( = 0.94892. 

One way of suppressing density fluctuations is via restricted random deposition (RRD): 
we divide the system into blocks, each consisting of 25 (or in some cases 26) sites, and 
depositing at random, within each block, a fixed number (24 or 25) of particles. In this way 
we generate a height distribution very similar to that of standard RD, but with virtually 
no density fluctuations on scales > 50 sites. Our conclusions from RRD simulations are 
as follows. First, the stationary activity density is identical, to within uncertainty, to that 
found previously. Thus we may continue to use f3/v± — 0.239 in the scaling analysis. Figure 
10, a scaling plot of the activity density (as in Fig. 8), shows a reasonable data collapse at 
long times, whereas the initial phase does not collapse. (It does collapse in an unrescaled 
plot, showing that the early stage is a size-independent, presumably purely local process.) 
Once again z has been chosen to optimize the collapse at long times, which in this case 
gives z = 1.47(6). The final approach to the stationary state appears to follow a power 
law, p a ~ t~ e , with 9 = 0.16(1). Analysis of the interface representation yields z = 1.57(3), 
a = 1.13(2) and ftw = 0.77(3), consistent with the scaling relation [3 W = a/z. (The growth 
of the survival time with L suggests yet another value, z ~ 1.3.) Thus for RRD, all of 
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the exponents describing transient properties are different from those found in RD studies. 
The exponent z has decreased to a value consistent with stationary relaxation, perhaps due 
to suppression of slowly relaxing, long- wavelength density variations. The sum 6 + j3 w , 
however, is now somewhat smaller than unity. 

A second kind of initial configuration with reduced density fluctuations may be termed 
maximum activity (MA). In this case all sites are either empty or have height 2 (except for a 
single site, in case N is odd), and the spacing between sites with z = 2 is fixed at 1 or 2 sites, 
in a regular pattern, so the average density is constant on a scale of 50 or so sites. While the 
density is again uniform on large scales, the number of active sites is as large as possible (as in 
usual initial-decay studies), and consequently the height distribution is quite different from 
that in the RD or RRD configurations. For MA IC's we find a = 1.12(4), (3 W = 0.78(1), and 
6 = 0.168(6). From the collapse of p a (t) and of W 2 (t), we obtain, respectively, z = 1.47(4) 
and z = 1.49(2), while the survival time data again yield z ~ 1.3. Thus our results for the 
two kinds of initial conditions, RRD and MA, with suppressed density fluctuations are in 
good agreement. This suggests that the precise form of the initial height distribution does 
not affect asymptotic scaling behavior. 

A general conclusion from these studies is that suppressing density fluctuations reduces 
the apparent value of the exponent z, suggesting that the unusually slow relaxation observed 
using random deposition initial conditions is due to long-wavelength density variations. The 
exponents a and (3w decrease as well, indicating that such density variations also provide 
a source of roughness for the interface dynamics. Further studies are required to determine 
whether the interface-growth exponents reported here are typical of an entire class of initial 
conditions (i.e., with reduced density fluctuations), or if their values may depend on other 
details of the initial preparation. Until such studies are performed, we cannot be sure 
that the exponents describing roughness and initial relaxation are well defined, describing 
asymptotic evolution, or merely represent apparent power laws over some limited range 
of times and system sizes. In none of the cases studied do we observe the simple initial 
relaxation of p a typical of the contact process. This may be due to the impossibility of 
preparing a completely uniform, high-activity initial state, and/or to the conservation of 
particles. 



E. Ergodicity 



In a recent study we found strong nonergodic effects in the BTW sandpile |T2|]. Given 
the randomness in the Manna sandpile rules, we do not expect such effects here. To test 
this hypothesis we plot (see Fig. 11) the active-site density p a versus time in four individual 
trials in a system of 1000 sites near ( c . The curves for p a in the various trials are completely 
intermixed; their ranges of fluctuation about a common mean are roughly equal. This is in 
complete contrast to the BTW case, where a similar plot shows distinct mean values, and 
very different ranges of variation of p a , in each trial |12| . Thus the Manna sandpile appears 
to be ergodic, in the sense that a time-average of p a in a single realization is independent 
of the initial condition (after a short transient), and yields the same value as the stationary 
average of p a over many trials. Nonergodicity in the BTW model has been associated with 
"hidden" conservation laws, i.e., a memory of the initial configuration that persists over the 
entire trial. The Manna sandpile seems to be free of such conservation laws: while we can 
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not rigorously exclude their presence, their effect must in any event be much weaker than 
in the BTW case. This is of course consistent with the fact that the BTW dynamics is 
deterministic, while the redistribution of particles is stochastic in the present case. (As a 
simple example of this, recall that a nearest-neighbor pair with both sites empty (z = 0) is 
forbidden in the stationary state of the BTW sandpile, whereas no such restriction exists 
for stochastic toppling rules.) 



IV. DISCUSSION 

We studied the scaling behavior of a one-dimensional fixed-energy sandpile with the same 
local dynamics as the Manna model. The model exhibits a continuous phase transition 
between an absorbing state and an active one at a critical particle density ( c = 0.94885(7). 
In the interface (integrated activity) picture this translates into a depinning transition with a 
fixed force (fixed-energy model) or a force that is ramped up at an infinitesimal rate (driven 
case) 0. 

The phase transition in the one- dimensional stochastic sandpile is characterized by the 
critical exponents (3 = 0.42(2) and v± = 1.81(1), which differ significantly from those asso- 
ciated with directed percolation (/? = 0.2765, u± = 1.0968) and linear-interface depinning 
((3 = 0.25(3), u± ~ 1.3). While absorbing-state phase transitions are expected to fall gener- 



ically in the directed percolation universality class |$7],[^|], it is reasonable to exempt the 
Manna model from this rule, due to local conservation of particles; this conservation law 
is expected to alter the universality class. In fact, studies of various models with the same 
local conservation law as the Manna sandpile, in dimensions d > 1, indicate a new, common 



universality class for models sharing this feature p2| . [21| . [22 l . 

Studying the interface representation of the model, we obtain the roughness exponent a = 
1.48(2) and growth exponent f3w = 0.86(2), which should be compared with a = 1.33(1), 
Pw = 0.839(1) for DP and a = 1.25(1), (3 W = 0.88(2) for LIM. Study of the height-height 
correlation function yields the local roughness exponent a/ oc = 0.59(3); the corresponding 
DP value is 0.63(3) pOfl . Changing the initial condition to suppress long- wavelength density 
fluctuations, we obtain a = 1.13(2) and j3w = 0.77(2), which again differ from the values 
associated with DP and LIM. Comparing these results with the apparent agreement of the 
Manna exponents with those of the LIM in two dimensions implies that either the rough 
numerical equivalence is fortuitous, or that the noise r in the interface equation has a 
fundamentally different structure depending on the dimension. It is worth remarking that 
the measured roughness exponent is rather close to a = 3/2, which is the value one expects 
if only the columnar component of the noise is relevant . The other exponents, however, 
seem to be far from the columnar disorder universality class (i.e., u±_ = 2, z = 2, /3 = 1, 

^ = 3/4) m. 



In linear interface models, translational invariance of the noise can be used to derive the 



scaling relation (2 — a)z/j_ = 1 |40|. This relation does not appear to be verified by our 



numerical results. Notice that in the context of interface depinning z can be linked to the 



other exponents by a scaling relation |^0[ which reads z = /3/v± + a. Inserting the values 
of /?, u± and a measured in simulations, we obtain z = 1.7, which is consistent with the 
exponent obtained analyzing the scaling of the width, for RD initial conditions. Our results 
for initial conditions with reduced fluctuations do not, however, follow this scaling relation. 
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Understanding this in the interface picture by considering the height-dependent quenched 
noise term is a challenge. 

Our results for the dynamic exponent z are conflicting. We have studied various defi- 
nitions of the relaxation or correlation time: that for the decay of the activity autocorre- 
lation in the stationary state (yielding z = 1.49(3)), one associated with the survival time 
(z = 1.42(1) for RD initial conditions, z ~ 1.3 for RRD and MA IC's), and those associated 
with the surface roughness and the decay of activity (z = 1.70(5) for RD IC's, z = 1.52(5) 
for initial configurations with suppressed density fluctuations). At a simple absorbing-state 
phase transition, all of these relaxation times are of course governed by the same exponent. 
Whether the unusually large z associated with time-dependent behavior using random de- 
position derives, as we have suggested, from slow relaxation of long-wavelength density 
fluctuations in the initial configuration, is a subject for further investigation. Our studies of 
modified initial states support this notion, but to reach a definitive conclusion other classes 
of initial conditions will need to be examined. Indeed, we have noted that, given their ap- 
parent dependence on the initial state, it is not clear whether the exponents z, a, /3w and 
9 are well defined, i.e., whether the asymptotic behavior follows simple power laws. 

We note that the present model does not exhibit the strong nonergodic effects observed in 
the fixed-energy version of the BTW sandpile. The relaxation of the mean height {zi(t)) from 
its initial value to the average, £, follows a power law with an exponent « 1/2. We find good 
evidence for finite-size scaling, in contrast with most driven sandpiles |34]-p6l . In summary, 
we have identified a one- dimensional sandpile model that exhibits an absorbing-state phase 
transition as the relevant temperature-like parameter (the energy density) is varied. It 
appears to be the "minimal model" for absorbing-state phase transitions belonging to a new 
universality class associated with a conserved density. Preliminary studies indicate that 



the driven version of the model exhibits scale-invariant avalanche statistics [0,0- We may 
therefore hope that analysis of the driven model, and of spreading of activity at the critical 
point of the fixed-energy system, will permit us to establish detailed connections between 
scale invariance under driving and the underlying absorbing- state phase transition. 

Let us finally stress that while in higher dimensions the linear interface universality class 
and that of systems with absorbing states in the presence of a conserved static field seem to 
coincide, the results of this paper show that this equivalence appears to be broken in d — 1. 
It will be interesting to study other one- dimensional systems with absorbing states and an 
order parameter coupled to a static conserved field P, p| , ^2| in order to compare the critical 



exponents and anomalies with those reported in this paper. 
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FIGURE CAPTIONS 



FIG. 1. Stationary active-site density versus energy density (. The points represent extrap- 
olations of data for L = 100 - 5000 to the L — > oo limit. 

FIG. 2. Stationary active-site density versus system size. From bottom to top, ( = 0.948, 
0.94857, 0.94864, 0.94874, 0.9488, 0.94892, 0.949, and 0.95. The inset shows the curvature 
b of the log-log plot as a function of (, for L > 1000. The straight line is a least-squares 
linear fit. 

FIG. 3. Inset: stationary autocorrelation function C(t) versus t at £ = 0.9488, for (left to 
right) L = 625, 1250, 2500, 5000, and 10 4 . In the main graph the data are plotted versus a 
rescaled time t* defined in the text. 

FIG. 4. Scaling plot of the active-site density versus A = ( — ( c . Symbols: +: L = 100; •: 
L = 200; x: L = 500; o: L = 1000; □: L = 2000; O: L = 5000. 

FIG. 5. Scaling plot of W 2 = W 2 /L 2a versus t = t/L z , for C = 0.9489, using a = 1.41, and 
z = 1.645(13). System sizes (top to bottom) L = 1253, 2506, 5012, 10024, and 20048. 

FIG. 6. Scaling plot of saturation width versus lifetime in individual trials, ( = 0.9490, 
L = 400-6400. 

FIG. 7. Saturation height versus lifetime in individual trials, ( = 0.9490, L = 400-6400. 

FIG. 8. Scaled active-site density p* = L l3 ^ u± p a versus scaled time t* = t/L z for ( = 0.9489, 
system sizes L = 1253,..., 20048 as indicated. 

FIG. 9. Mean height (z(t)\z(0)) of sites with initial height z(0), in at system of 1400 sites 
at ( c , averaged over 2000 trials. From bottom to top: z(0) = 0, 1, 2, 3, and 4. The inset is 
a plot of In \ (z(t)\z(0)) — (\ versus Int. 

FIG. 10. Scaled active-site density versus scaled time for ( = 0.9489 as in Fig. 8, but for 
RRD initial conditions. 

FIG. 11. Active-site density p a (t) in four different trials at ( — 0.95 in a system of 1000 
sites. 
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